//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Quadrature.hh
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Quadrature class definition.
 * \note   Copyright (C) 2011 Jeremy Roberts.  Derived from Denovo.
 */
//---------------------------------------------------------------------------//
// $Rev:: 123                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-08-25 12:30:43 +0000 (Thu, 25 Aug 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#ifndef QUADRATURE_HH
#define QUADRATURE_HH

#include <iostream>
#include <vector>

#include "angle/AngleTraits.hh"
#include "tran/Definitions.hh"
#include "tran/Traits.hh"
#include "../utilities/DBC.hh"

namespace slabtran
{

//===========================================================================//
/*!
 * \class Quadrature
 * \brief A class to encapsulate the angular discretization.
 *
 * The Quadrature class provides services related to the angular
 * discretization scheme.  It creates a set of quadrature directions
 * (abscissas) and weights associated with a particular quadrature scheme
 * specified by the calling routine.
 *
 * The angles in the quadrature are defined:
 * \f[
 \mathbf{\Omega} = \Omega_x\mathbf{i} + \Omega_y\mathbf{j} +
 \Omega_z\mathbf{k}\:.
 * \f]
 * Where the components of \f$\mathbf{\Omega}\f$ are defined
 * \f[
 \Omega_x = \sin\theta\cos\varphi\:,\
 \Omega_y = \sin\theta\sin\varphi\:,\
 \Omega_z = \cos\theta\:.
 * \f]
 * Here, \f$\theta\f$ is the polar-angle and \f$\varphi\f$ is the azimuthal
 * angle. Following common convention, we use the symbols \f$(\mu,\eta,\xi)\f$
 * for the components of \f$\mathbf{\Omega}\f$:
 * \f[
 \mathbf{\Omega} = \mu\mathbf{i} + \eta\mathbf{j} + \xi\mathbf{k}\:.
 * \f]
 * The solid angle is integrated over a unit-sphere:
 * \f[
 \int_{4\pi}d\mathbf{\Omega} =
 \int_{0}^{2\pi}\int_{0}^{\pi}\sin\theta\ d\theta d\varphi = 4\pi\:.
 * \f]
 * The first moment of the integral is:
 * \f[
 \int_{4\pi}\mathbf{\Omega}\ d\mathbf{\Omega} = 0\:.
 * \f]
 * The second moment of the integral is a tensor:
 * \f[
 \int_{4\pi}\mathbf{\Omega\Omega}\ d\mathbf{\Omega} =
 \left(\begin{array}{lll}
 \frac{4\pi}{3} & 0 & 0 \\
   0 & \frac{4\pi}{3} & 0 \\
   0 & 0 & \frac{4\pi}{3} 
 \end{array}\right)\:.
 * \f]
 *
 * Quadrature stores the following data:
 * \f[
 (\mu_a, \eta_a, \xi_a: w_a)\ a=1,n\:,
 * \f]
 * where \f$ n\f$ is the number of angles in the set and \f$ w_a\f$ is the
 * weight of each angle.  Using the quadrature the integrals shown above are
 * expressed:
 * \f[
 \int_{4\pi}d\mathbf{\Omega} = \sum_{a=1}^{n}w_a = C\:.
 * \f]
 * The user is allowed to choose the normalization constant, \f$ C\f$.  We use
 * \f$ C=4\pi\f$ in denovo to be consistent with the analytic integrals.
 *
 * \note Many Nuclear Engineering texts use a left-handed coordinate-system in
 * which \f$\mu=\cos\theta\f$.  This allows 1-D slab problems to be defined in
 * \f$ x\f$ instead of \f$ z\f$.  We choose to use a more
 * "mathematically-standard" coordinate system in which the polar angle
 * (\f$\theta\f$) is defined from the \f$ z\f$-axis.  Also, it is common in
 * many codes to set the normalization of the quadratures to 1 as this removes
 * the annoying \f$ 1/4\pi\f$ factors everywhere.  We choose to let the
 * quadrature integrate to \f$ 4\pi\f$, thus being consistent with the
 * analytic integral of \f$ d\Omega\f$.
 */
//===========================================================================//

template <class D>
class Quadrature
{

public:

  // Typedefs.
  typedef std::vector<double>               Vec_Dbl;
  typedef typename QuadTraits<D>::Angles    Angles;

  // Constructor.
  Quadrature(int sn_order, double norm);

  // Virtual destructor.
  virtual ~Quadrature() = 0;

  // >>> VIRTUAL INTERFACE

  //! Return the \f$ S_N\f$ order.
  int sn_order() const
  {
    return d_sn_order;
  }

  // >>> ACCESSORS

  //! Return the angles.
  const Vec_Dbl& angle(int angle = 1) const;

  //! Return the angle for a given ordinate.
  double angle(int angle, int m) const;

  //! Return the dimensionality.
  int dimensionality() const;

  //! Get the weights.
  const Vec_Dbl& weight() const
  {
    return d_wt;
  }

  //! Return the weight for a given ordinate.
  double weight(int m) const
  {
    Require (m >= 0 && m < d_wt.size());
    return d_wt[m];
  }

  //! Return the number of angles.
  int number_angles() const
  {
    return d_wt.size();
  }

  //! Return the normalization.
  double norm() const
  {
    return d_norm;
  }

  // Display the quadrature.
  virtual void display() const;

  // Display the quadrature name.
  virtual void name() const = 0;// { std::cout << "BASEQUAD..." << std::endl; } ;

protected:

  // >>> DATA

  // Sn-order.  Can be ambiguous!
  const int d_sn_order;

  // Normalization.
  double d_norm;

  /// Quadrature directional cosines.
  Angles d_angles;

  /// Quadrature weights.
  Vec_Dbl d_wt;

};

} // end namespace slabtran

#endif // QUADRATURE_HH
//---------------------------------------------------------------------------//
//              end of Quadrature.hh
//---------------------------------------------------------------------------//
